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O ABSTRACT 
(N 

<U I We develop an unconditionally stable numerical method for solving the coupling 

between two fluids (frictional forces/heatings, ionization, and recombination), and in- 
vestigate the dynamical condensation process of thermally unstable gas that is pro- 
' r-| ' \ vided by the shock waves in a weakly ionized and magnetized interstellar medium by 

using two-dimensional two-fluid magnetohydrodynamical simulations. If we neglect 
O I the effect of magnetic field, it is known that condensation driven by thermal instability 

^ I can generate high density clouds whose physical condition corresponds to molecular 

clouds (precursor of molecular clouds). In this paper, we study the effect of mag- 
netic field on the evolution of supersonic converging HI flows and focus on the case 
>• ' in which the orientation of magnetic field to converging flows is orthogonal. We show 

\o ■ 

00 

'nI" ■ density and high column density clouds, but instead generates fragmented, filamentary 

HI clouds. With this restricted geometry, magnetic field drastically diminishes the op- 
(3 ■ portunity of fast molecular cloud formation directly from the warm neutral medium, 

^ ■ in contrast to the case without magnetic field. 

> 

. 1. Introduction 

It is widely known that radiative cooling and heating in the interstellar medium (ISM) leads 
to a stable cold dense medium and a warm diffuse medium (Field, Goldsmith & Harbing 1969; 
Wolfire et al. 1995, 2003). These stable phases are called the cold neutral medium (CNM), and the 
warm neutral medium (WNM). The CNM is observed as HI clouds {n ~ 10^"^ cm~^, T ~ 10^ K), 
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and cold neutral clouds that contain substantial amounts of H2 are called molecular clouds (n ~ 
10^~^ cm~^, T ~ 10 K). Since molecular clouds are the sites of star formation, the knowledge of 
how molecular clouds are formed is crucial to understanding the initial conditions of star formation. 

Thermal instability (TI) is the most promising formation mechanism of the CNM. Many au- 
thors studied the dynamical condensation process of the ISM driven by TI. Koyama & Inutsuka 
(2000, 2002), and Inutsuka, Koyama, & Inoue (2005) studied the CNM formation as a result of 
TI in the layer compressed by shock propagation. Hennebelle & Perault (1999), Audit & Hen- 
nebelle (2005), Hennebelle & Audit (2007), Heitsch et al. (2005, 2006), and Vazquez-Semadeni 
et al. (2006, 2007) studied an analogous process in supersonic converging flows, which essentially 
include two shocks. These studies showed that TI can generate molecular clouds, by piling up the 
WNM. 

However, the above-mentioned studies do not include the effect of magnetic field. The dy- 
namical formation process of the CNM with the effect of magnetic field is studied by Hennebelle & 
Perault (2000) using one-dimensional, ideal MHD simulations. They showed that magnetic pres- 
sure prevents TI, if the initial angle between the magnetic field and the fluid velocity is larger than 
20°-30° with a few microgauss initial field strength. Recently, Inoue, Inutsuka, & Koyama (2007) 
have shown that if the effect of ambipolar diffusion is taken into account, TI can create the CNM 
even with a microgauss magnetic field which is orthogonal to the orientation of the condensation. 

In Inoue, Inutsuka, & Koyama (2007), we used a thermally unstable equilibrium as an initial 
condition and performed one-dimensional numerical simulations. However, thermally unstable 
gas is likely to be supplied by the shock compression of the ISM, and in multiple dimensions, 
growth of TI along the magnetic field seems to be faster than the mode that is perpendicular to the 
magnetic field. In this paper, to investigate how the CNM is formed under more realistic situations, 
we study the dynamical condensation process of thermally unstable gas that is supplied by the 
shock compression of the WNM employing two-dimensional two-fluid magnetohydrodynamical 
simulations. 

This paper is organized as follows. In §2 we provide the basic equations of the weakly ionized 
interstellar medium and assumptions used in this paper. In §3 we describe a method for solving 
the basic equations numerically, in which an unconditionally stable scheme that solves the cou- 
pling between two fluids is proposed and tested. The results of two-fluid simulations of TIs in 
weakly ionized media are shown in §4. Finally, in §5 we summarize our results and discuss their 
implications. 
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2. The Model 

2.1. Basic Equations 

We consider the equations of two-dimensional hydrodynamics (HD) and magnetohydrody- 
namics (MHD) that include the effects of ionization, recombination, and coUisional coupling be- 
tween neutral and ionized gases (Draine 1986). The effects of radiative cooling, heating, and 
thermal conduction are taken into account in the neutral fluid part. 
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where e = p/(7 — 1) + pv'^/2 is the total energy per volume and V = {d/dx, d/dy, 0) is the 
two-dimensional spatial derivative operator. The subscripts n and i denote the neutral and charged 
components, respectively. In the neutral component, we take into account hydrogen and helium 
atoms. In the charged component, we take into account hydrogen, helium and carbon ions. The 
ionized hydrogens and heliums are described with the subscript p, and the ionized carbons are 
described with the subscript CII (pi = Pp + pen)- We assume that the density of the ionized carbon 
is proportional to the density of the neutral gas pci = AcPn'mcn/fTT-n^ where we use = 3x 10"'^. 
In addition, we impose the ideal gas equation of state pn = k-^pnT / rrin, and we use the isothermal 
approximation between the neutral gas and the ionized gas p^ = pn nin Pi/m-J pn, where the mean 
charged particle mass is determined by mi = (pp + pci) / {up + ncn) • We use the specific heat ratio 
7 = 5/3 and the mean particle masses rrin = rrip = 1.27mpi.oton (i-C-^ we assume for convenience 
that 91% of neutrals and light ions are, respectively, hydrogens and hydrogen ions, and that the 
remaining 9% of neutrals and light ions are, respectively, heliums and helium ions). The drag 
coefficients Ap and Aqji, the recombination coefficient a, the ionization coefficient ^, and the 
thermal conductivity k are chosen to simulate the typical ISM whose properties are listed in Table 
1. 
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2.2. Cooling Function and Its Characteristics 



We take the following simplified net cooling function: 
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This is obtained by fitting to the various heating (F) and cooling (A) processes considered by 
Koyama & Inutsuka (2000). The resulting thermal-equilibrium state defined by C{n^, T) = is 
shown in Figure [T] as a thick solid line. We also show the contour of the timescale of the cooling 
and heating tcooi = ksT / {m^C (7 — 1)} as broken lines. 

Isochorically cooling gas and isobarically contracting gas, which can be generated by shock 
compression, are thermally unstable, if these gases satisfy the Balbus criterion [d{C/T)/dT]p < 
(Balbus 19950 In Figure HI we show the unstable region as a dotted region. 



We use the operator- splitting technique for solving the basic equations which are split into 
four parts: (1) the ideal HD part for neutral gas, (2) the ideal MHD part for ionized gas, (3) the 
cooling, heating and thermal conduction part, and (4) the ionization, recombination, and frictional 
forces/heating part. 

Methods that solve parts 1, 2, and 3 are well established. We solve the ideal HD part by 
employing the second-order Godunov method (van Leer 1979). The ideal MHD parts are solved by 
employing the second-order method developed by Sano, Inutsuka & Miyama (1999), in which the 
HD parts with the magnetic pressure term are solved using second-order Godunov method, and the 
magnetic tension term and the induction equation are solved using the method of characteristics for 
Alfven waves (Stone & Norman 1992) with the constrained transport algorithm (Evans & Hawley 
1988). We use the explicit time integration for the cooling, heating and thermal conduction, which 
has second-order accuracy in time and space. 



3. Numerics 



3.1. Numerical Schemes 



'The linear stability analysis of isochorically cooling gas was done by Schwarz, McCray, & Stein (1972) and that 
of isobarically contracting gas was done in Appendix B of Koyama 8c Inustuka (2000). 
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In the parts 1 and 2, equations are solved in conservative fashion. Thus, the total mass, 
momentum, and energy of neutral gas and the total mass and momentum of ionized gas conserve 
exactly in the absence of source terms. The internal energy of ions is determined by equating 
temperatures of neutral and ionized gases. This approximation is reliable when the ionization 
degree is small and the relaxation time is short. In HI regions considered in this paper, both 
conditions are fulfilled. It is known that, when the plasma /5 is small (< 0.01), the conservative 
MHD scheme that uses total energy (E = e + B"^ / Svr) as a basic variable to update cannot describe 
thermal energy accurately, because thermal energy should be calculated by subtracting magnetic 
energy (fi^/Svr) and kinetic energy {pv"^ /2) from the total energy {E). In the simulation of weakly 
ionized plasma employing a two-fluid conservative method, we would encounter this problem even 
when the (3 of the whole medium is not small, since the plasma (3 of ionized gas is always lower by 
a factor of the ionization degree. However, if we assume the temperatures of neutral and ionized 
gases are the same, we can perform the simulation without such difficulty, as in this paper. 

In part 4, we use the method presented in the next section which utilizes exact solutions of 
split equations and is unconditionally stable. Thus, part 4 does not restrict the time step, and the 
time step in our code as a whole is determined so that the following criteria are satisfied every- 
where: the CFL condition for the sound wave of neutral gas (Atn < 0.5 Ax/(cs + l^'nl)). the 
CFL condition for the fast wave of ionized gas (Ati < 0.5 A x/(cfast + the condition for 

avoidance over-cooling/heating (Atcooi < 0.2pn/[pn l-^^l (7 — 1) ]), and the stability condition for 
thermal conduction (Atcond < 0.5 (Ax)^ A;b/[ (7 — 1) n])- 



3.2. Source Term Solver 

3.2.1. Exact Maps 

In this section we present a scheme for solving the source part in the basic equations ([I])-®. 
Let us consider differential equations 

= 9iU)+92U). (10) 

^ = 92{f), (11) 
df 

= 9i{f), (12) 

where g is an operator which can include the spatial derivatives of /. If we introduce a map 
G{At, f) that gives a solution of a differential equation at t + At with second-order accuracy in 
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time, 

f{t + At) = G'i2(At,/(t)) + 0(At3), (13) 

f{t + At) = Gi(At,/(t)) + 0(At3), (14) 

fit + At) = G2(At,/(t)) + 0(At3), (15) 

then it is known that the operator Gu can be approximated by the combination of Gi and G2, 

Gi2(At, /) = Gf Gf Gf / + 0{At'). (16) 

where we have introduced a convention G{Ati, H{At2, /)) = G^^^ H^^'^ f. In the case of basic 
equations Q-®, the operator gi corresponds to HD, MHD, cooling, heating, and thermal con- 
duction operators, which are also solved by using further operator splitting explained in §3.1 , and 
the function g2 corresponds to other source terms. The approximation from equation (fT6l) remains 
valid, as long as the map G2 is more accurate than the second-order. In our case, the split source 
equations are the following. 
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We can derive the following exact map (t ^ t + At): 

Unit + At) = ran(t) exp~«^*, (24) 

hi{t + At) = +nn(t) (1 -exp-«^*), (25) 

- (+ j_ _ n^jt) + aAtnp{t) {nnjt) + np{t)} 

+ ^ " l + aAtnp(t) ' ^ 

Ut + At) = , (27) 
1 + a Atnp{t) 

Vr.it + At) = -^—^{v,(t)A,{l-exp-^^^+^-^^') 

(^2 + A, exp-(^^+^^)^*)}, (28) 
i^i(t + At) = ^^{i;„(t)A2(l-exp-(^^+^^)^*) 

+6,(«)(Ai + /l2exp-<^'+-*»>^')}, (29) 
p„(t + AO = pW + ~ ''^;ff {1 - exp--'-^-^^-'^'}. (30) 

where Ai = Ap pp + Acn Pci, ^2 = (^p Pp + ^ci Pen) Pn/pi, and A3 = pn {(^p rrip pp)/ (mn + 
rrip) + (v4cn "^ci Pen) / ("^n + "^ci)}: and we assume that the coefficient (a function of temper- 
ature) is constant during time step At. Even in the case where we do not use the equal temperature 
approximation between neutral and ionized gases, we can find the exact maps for pn and pi that 
were used in Inoue, Inutsuka & Koyama (2007). 

Since variables are evolved according to the exact solutions, this part is unconditionally stable 
and does not limit the time step. Furthermore, if G2 is a exact map, the time evolution {t 
t + S"^]^ At) can be expressed as 



£ Ij^ _|_ \ ^ At ) ^At„ Q^tn/^, ^At„_l/2 ^At„_i ^At„_l/2 

>> V / J 2 12 2 12 



1=1 



^Ati/2 ^M, ^^^^ ^ ^(AtS) (31) 

^At„/2 ^At„ ^(Ai„+At„_i)/2 ^At„_i ^(At„_i+At„_2)/2 
'^2 '^1 '^2 '^l '^2 

^iAt2+At,)/2 ^Ati ^Ati/2 ^^^^ ^ 0{At^). (32) 



In equation (I3T|) . there are 3 n operators, while they become 2n + 1 operators in equation (132 
Thus, by using the exact map, we can save computation time without loss of accuracy. 
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3.2.2. Tests for Alfven Waves with Plasma Drift 

In order to confirm the capability of our method, we perform test calculations with various 
collision frequencies. We also verify our method by a convergence test. The basic equations used 
in these tests are 

Qv 

Pn = A pn Pi {Vi,y - Vn,y), (33) 

Pi = Bx ^ + ^ Pn Pi {Vr,,y - Vi^y), (34) 

dBy _ dV.^y 

~ B^ ~n — • (33) 

Ot OX 

If we assume pn, Pn, and A are constants, these equations become linear differential equations 
that describe the evolution of Alfven waves in a uniform partially ionized medium. Assuming the 
solutions to be oc exp x — ujt), the characteristic equation can be written as (Kulsrud & Pearce 
1969) 



UO [UJ 



2 k'')+iv {{l + e)uo'' -clk'^e] (36) 



where u = Ap^h the collision frequency of a charged particle in a sea of neutrals, e = pj Pn is 

1 /2 

the ionization degree, and ca = B^/ p-J is the Alfven velocity in the weak-coupling limit. This 
characteristic equation contains three modes. In Figure [H we plot the real (top) and imaginary 
(bottom) parts of uj against the collision frequency u. The other parameters are chosen such that 
Bx = 1.0, Pn = 1.0 X 10^, Pi = 1.0, and k = 1.0, which are used throughout in test calculations. 
The solid lines are the branches of the Alfven waves that propagate in the weak (0 < < 1.96) and 
strong (u > 5.03) coupling regimes. The dashed lines show the other mode, i.e., a purely damping 
mode. In the weak-coupling limit (u 0), the phase velocity of the Alfven wave is given by 
Bx/pI^^ = 1-0, and in the strong-coupling limit (z/ oo), it is given by B,j,/{pi + Pn)^^^ = 
9.95 X 10-2. 

To demonstrate the propagation of Alfven waves, we set the initial condition by eigenfunc- 
tions of three modes, 

By{x) = Re [e^^] , (37) 
-Re[ujBy{x)], (38) 



Vn,y{x) = Vi^yix) + Re 



^{uj'-l)By{x) 



(39) 



We take the cu that is the solution of the characteristic equation (|36l) and is on the branch of the 
thick solid line in Figure [IJ In the test code, we use an operator- splitting technique based on 
equation (|32|) in which we use the MOC scheme (Stone & Norman 1992) and the exact map of the 
friction part based on equations (l28l) and (|29l) . The former one solves the propagation of Alfven 
waves in the absence of the friction term, and the latter one solves the friction force term. We use 
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a periodic domain that covers the x G [0.0, 2 vr] region. The time step At is determined from the 
CFL condition for the Alfven wave (in this case, At = vr/m, where m is the cell number used in a 
run). 

The second test is a convergence test for the case v = 10.0. Calculations are performed from 
t = 0.0 to 20.0, while the Alfven waves propagate about 0.28 wavelengths and are damped by a 
factor 0.37. Since basic equations are linear differential equations, we can make direct comparison 
between numerical and exact solutions. In Figure [H we show both the numerical solution with 
m = 128 and the exact solution at t = 20.0. The numerical solution agrees well with the exact 
solution. In Figure |3] we show the plot of average deviation from the exact solution att = 20.0, 

^^"^) = m ^ '^^'^ ~ ^^(^Oli=2o.o ' (40) 
1=1 

as a function of cell number m used in calculations, where By i is the magnetic field at the i-th cell 
center calculated numerically, and By{x) is the exact solution. This figure clearly shows that the 
second-order accuracy in time and space is achieved in our scheme based on the exact map and the 
transformation corresponding to equation (|32l) . 

In order to confirm the capability of our method from the weak-coupling regime to the strong- 
coupling regime, we examine the propagation of Alfven waves in a partially ionized medium 
with various collision frequencies. The propagation speed of the numerical wave can be de- 
termined by tracing the position where By is maximum. We measure the propagation speed 
and calculate numerical Re[uj] by using the time when the position of maximum By crosses 
X = 271, except the cases of 1.96 < i' < 5.03 in which the propagation speed is measured at 
the time t = 20.0. The damping rate (numerical Im[ct;]) is measured by using the time when 
Jq'^ \By(x,t)\dx/ Jq'' I 0)1 da: becomes smaller than 1/e. In Figuredlwe plot the numerically 
obtained cuhy m = 128 calculations as points. This figure clearly shows that our method can solve 
the two-fluid system in all the cases from the weak-coupling regime to strong-coupling regime. 



3.3. Initial Settings 

In order to study the generation of clouds as a result of TI, we consider the situation in which 
supersonic WNM flows collide with supersonic velocity, which will supply thermally unstable 
gas in the shocked slab. A rectangle computational domain is used whose side lengths are 30.0 
pc for the a:-direction (x = —15.0 to 15.0 pc) and 2.0 pc for the y-direction (y = —1.0 to 1.0 
pc). Two opposing WNM flows with density fluctuations initially collide head-on at the x = 
0.0 interface. Periodic boundary condition is imposed on boundaries at y = —1.0 and 1.0 pc, 
and WNM flows without fluctuations are constantly set at x = ±15.0 pc boundaries. Note that 



the density fluctuations are added only initially, and their properties are described below. We 
use nonuniform mesh. The region of a; = —1.0 to 1.0 and y = —1.0 to 1.0 pc is divided into 
2048 X 2048 cells, and the regions of x = -15.0 (15.0) to -1.0 (1.0) pc and y = -1.0 to 1.0 pc 
are divided into 512 x 2048 cells. Thus, the spatial resolutions are Ax = 9.8 x 10~^ pc in the 
higher resolution region. Ax = 2.7 x 10"^ pc in the lower resolution region, and Ay = 9.8 x 10"^ 
pc. In the higher resolution region, the Field condition (Koyama & Inutsuka 2004) is satisfied. 

We choose the density and pressure of the unperturbed WNM as p„ = 0.57 m,! and = 
3.5 X 10'^ /cb, which are the thermal equilibrium values of the cooling function (eq. dll). In order 
to excite TI, we add fluctuations in the initial density of the WNM as 



logp„,ini = I 1 + X] X] 



^ ^ ^'271 kyy 



''X 



sm 



Ly 



e{iy 



logpw, (41) 



where p„,cq is the neutral gas density of the unstable equilibrium, 6{l) is a random phase, and we 
use A = 0.1. By doing so, fluctuations whose scales are around A = 1 pc (the most unstable scale 
of TI) are generated in the shocked layer as a result of the WNM flows collision. The minimum and 
maximum density of the perturbed initial condition are Pn,min = 0.37 mn and Pn.max = 0.93 mn. 
The initial weakly ionized gas is set as an ionization equilibrium. 

The initial magnetic field strength is 2.0 /iG, which might be a typical value in the WNM. The 
initial orientation of the magnetic field is perpendicular to the initial inflow velocity, so that the 
effect of magnetic field becomes the strongest. The effect of the orientation of the magnetic field 
is discussed in §4. In order to contrast the effect of magnetic field, we also performed the same 
calculation except in the absence of magnetic field. 

We set the initial converging velocity at 20.0 km s^^, with which the Mach number is M = 
I'^^xl/cfast = 2.1 for the unperturbed gas of the magnetized case and M = \vx\/cs = 2.4 for that 
of the unmagnetized case. We also performed several calculations decreasing the initial velocity 
down to M = 1 both with and without magnetic field, and confirmed that the results shown below 
were not qualitatively changed. Note that as shown in Hennebelle & Perault (1999), if the speeds 
of the converging flows are subsonic, thermal condensation due to TI does not take place. 



4. Results 

The evolutions of the shocked WNM are very different in magnetized and unmagnetized 
cases: 



Unmagnetized case. 



- 11 - 



At first, as indicated in Koyama & Inutsuka (2000), shock-compressed gas gradually contract 
almost isobarically, until the gas satisfies the Balbus criterion. Once the gas becomes thermally 
unstable, overdensity perturbations begin to condense exponentially in time, until these regions 
reach a stable CNM phase. In Figured we show the evolutionary track of the highest density gas 
in the computational domain as a red line until it reaches the stable CNM phase, which is at t = 0.5 
Myr. 

The density distribution at time t = 2.43 Myr is shown in the left panel of Figure [51 In- 
side the shocked slab, CNM cloudlets are generated, and the interwoven medium composed of 
the cloudlets, unstable gas, and WNM is formed. In the core of the generated cloudlets, number 
density reaches n > 1 x 10'^ cm"'^. Since, these small cloudlets move randomly with supersonic ve- 
locity dispersion, as first pointed out by Koyama & Inutsuka (2002), the shocked slab has complex 
geometry (see also Heitsch et al. 2005). 

Despite the fact that the shocked gas is cooled until it reaches the CNM phase, the thermal 
pressure of the CNM is maintained at a high value (p/ ^ 1 x 10^ K cm^'^). This is because the 
CNM and surrounding thermally unstable gas slab are bound by the ram pressure of the converging 
flows. When the cooling rate in the slab balances the thermalization rate of the shocks, outward 
propagating shocks become quasi-standing shocks. This transition takes place at t ~ 0.5 Myr in 
this calculation, which is roughly simultaneous with the time of the CNM formation. Thus, as long 
as the converging flows continue, high-density and high-pressure conditions in the CNM are held, 
and the mass of the CNM continues to increase owing to the accretion of the "fresh" unstable gas 
which is continuously supplied by the quasi-standing shocks. Such a CNM would evolve into a 
molecular cloud. Similar results are also reported in the papers introduced in §1. 

• Magnetized case. 

At first, shock-compressed gas isochorically cools, decreasing its thermal pressure, until the 
gas satisfies the Balbus criterion. Once the gas becomes thermally unstable, overdensity perturba- 
tions begin to condense. During the growth of the TI (condensation), the thermal pressure contin- 
ues to decrease due to the cooling. Since the growth timescale of the most unstable scale of TI is 
equal to the timescale of the cooling (Field 1965), by the time condensing regions reach a stable 
CNM phase, the thermal pressure has fallen p/k^ — 2 x 10'^ K cm^^, which is on the order of the 
average pressure of the ISM. In Figure|4]we show the evolutionary track of the highest density gas 
as a blue line until it reaches the stable CNM phase, which is at t = 1.6 Myr. This evolution is 
quite different from the case of the unmagnetized result. This is because magnetic pressure in the 
shocked slab can sustain the ram pressure of the converging flows, and the thermal pressure does 
not have to be high in order to oppose the ram pressure. Thus, in contrast to the unmagnetized case. 
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the shocks never become standing shocks and the width of the post shock layer increases almost 
linearly with time. 

The density distribution at time t = 2.43 Myr is shown in the right panel of Figure[5l Because 
the decrease of the thermal pressure leads to longer cooling times (see contours of cooling time 
in Figure HJ, the nonlinear growth of the TI in the magnetized case needs a longer time than the 
unmagnetized case. At this time, most of dense clumps have reached the stable CNM phase in the 
region shown in the figure. We also show the density distribution (Fig. 6, top) and magnetic field 
lines (Fig. 6, bottom) of the whole simulation domain in Figure [6] at the same time depicted in the 
right panel of Figure [51 which indicates that the propagation of shocks supply thermally unstable 
gas to an extended region. It seems that oblong grid cells used in lower resolution regions do not 
show any artificial effect. In the post-shock region (region where > 1.0 cm~^), the average 
magnetic field strength is \B\ = 11.7 fxG. 

The density of the CNM clumps is typically n ~ 3 x 10^ cm"'', which corresponds to that 
of HI clouds. The shapes of the CNM clumps are filamentary, because the growth of TI due 
to the motion along the magnetic field is faster than the motion perpendicular to the magnetic 
field. We also performed a similar calculation for the ideal MHD equations subject to the same 
cooling function and thermal conduction with the same initial and boundary conditions. The result 
was very similar to the two-fluid case. In a one-dimensional case, as shown in Inoue, Inutsuka 
& Koyama (2007), ambipolar diffusion is effective for the scale of the TI with the presence of 
microgauss magnetic field. In a two-dimensional case, however, the mode along the magnetic field 
dominates the nonlinear growth, and the role of ambipolar diffusion becomes less important at least 
in the formation regime of the CN]v|^. Therefore, even if the converging flows last long, TI cannot 
pile up material in the direction parallel to the converging flows, and it is difficult to generate a 
high-density cloud that resembles molecular clouds, since the propagating shocks supply thermally 
unstable gas to an extended region. 



5. Summary and Discussion 

We have developed an unconditionally stable numerical method for solving the coupling be- 
tween two fluids, which can adequately solve the evolution of a partially ionized medium from 
weak-coupling to strong coupling regimes. By using two-dimensional two-fluid MHD simulations 
based on this method, we have investigated the dynamical formation process of the CNM in both 
the magnetized and unmagnetized media. Our findings are as follows. 



^The role of ambipolar diffusion may be important through frictional heating in the regime of late time evolution 
of the two-phase medium as studied in Hennebelle & Inutsuka (2006) 
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In the unmagnetized case, the converging flows of the WNM create a complex slab bound by 
quasi- standing shocks in which turbulent cloudlets of the CNM are generated as a consequence of 
TI. Because the thermal pressure in the slab has to sustain the ram pressure of the flows, the CNM 
whose pressure is higher than the average pressure of the ISM and whose density is n ~ 10^ cm""^, 
which correspond to the observational characteristics of molecular clouds, is inevitably formed. 
Thus, as long as the converging flows continue, the column density of the CNM increases, and a 
molecular cloud would be generated in the slab. 

In the magnetized case, magnetic pressure inside the shocked region can sustain the ram 
pressure of the flows. This leads to two important effects on the evolution of the shocked region. 
One is that, because the thermal pressure does not have to be high in order to oppose the ram 
pressure, it decreases due to cooling. This makes the physical conditions of the resulting CNM 
similar to those of HI clouds. Another one is that the shocks continue to propagate, supplying 
thermally unstable gas to an extended region. This makes it difficult to generate high-column 
density clouds in a compact region like molecular clouds. This is because ambipolar diffusion is 
not effective globally, and TI cannot pile up material in the direction parallel to the converging 
flows. 

In this paper, as a first step, we examined the situation in which the initial angle between 
the magnetic field and the velocity of the colliding flows is 90°, which maximizes the effect of 
magnetic field. If the initial angle is small enough, it seems possible to generate the CNM whose 
physical condition corresponds to molecular clouds, as in the unmagnetized case. How large is 
the critical angle that changes the character of the resulting CNM? In a one-dimensional case, 
according to Hennebelle & Perault (2000), magnetic pressure prevents condensation, if the initial 
angle is larger than 20°-30° with a few microgauss initial field strength (see Fig. 10 of their 
paper). In more realistic two-dimensional situation, as we have shown in this paper, the mode 
of TI along the magnetic field always generates the CNM that corresponds to a HI cloud. Thus, 
their conclusion needs to be modified. However, their result tells that if the initial angle is smaller 
than their critical angle, magnetic pressure cannot be strong enough to support the ram pressure. In 
such a case, the resulting shocked slab would have a physical condition that can generate molecular 
clouds (precursor of molecular clouds), provided the flows last long. 

However, if this conjecture is right, smallness of the critical angle will drastically diminish 
the opportunity of fast molecular cloud formation directly from the WNM whose time scale is on 
the order of growth time scale of TI, or cooling time, (~ 1-10 Myr). Thus, there is a possibility 
that most of the molecular clouds are formed from HI clouds, as suggested from observations of 
clouds in Local Group galaxies (Fukui 2007; Blitz et al. 2007), via more complex processes (e.g., 
experiences of TI more than once due to successive shocks of supemovae and possibly through 
the effect of potential well of stellar spiral waves or self-gravity) rather than the direct pile-up of 
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WNM as a consequence of TI. The diffuculty of the fast formation of molecular clouds in the 
presence of magnetic field perpendicular to the inflow is also discussed from another point of view 
in Hartmann et al. (2001). 

In order to confirm this possibility, the following studies are necessary: (1) similar calculations 
of converging flows varying the angles between converging flows and magnetic field to show the 
smallness of the critical angle even in multiple dimensions; and (2) a global study of the magnetized 
ISM to learn about the typical angle between the mean magnetic field and converging flows or 
interstellar shocks. In subsequent papers, we will present studies relevant to these points. 

Note that the discussion shown in this section is based on one-/two-dimensional simulations. 
However, especially in MHD, dynamical behavior in 3D can be different from that in 2D (see, e.g.. 
Stone & Gardiner 2007 for Rayleigh-Taylor instability in a magnetized medium). Thus, we should 
bear in mind the possibility that the difference due to the dimensional restriction might change the 
outcome. 
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Table 1. Sources 



Coefficient 


Note 


Ref. 




3.4 X 10^^ (T/8000 K)°-^^^ cm^g-is^^ 


1 




H"*'+H collisions 






8.5 X 10^2 cm^g-^s"^ 


2 




Cn+H collisions 




K 


2.5 X lO^T^'S erg cm-^s-^K-i 


3 




H+H collisions 




a 


H+ and He+ recombinations 


4 


e 


1.0 X 10-16 s-i 


5 




H and He ionizations 





References. — (1) Glassgold et al. 2005; (2) Draine te al. 
1983; (3) Parker 1953; (4) Shapiro & Kang 1995; (5) Wolfire 
et al. 1995. 
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coUision frequency • v = A 



Fig. 1 . — Real (top) and imaginary (bottom) parts of lu against the collision frequency u in the case 
Bx = 1.0, pn = 1-0 X 10^, Pi = 1.0, and k = 1.0. The solid lines are the branches of the Alfven 
waves that propagate in the weak (0 < < 1.96) and strong (v > 5.03) coupling regimes. The 
dashed lines show the other mode, i.e., a purely damping mode. Propagation speeds (numerical 
Re[c<j]) and damping rates (numerical Im[ci;]) obtained from test runs with cell number m = 128 
are plotted as points in the cases that u = 0.5, 1.0, 1.5, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, and 
10.0. 
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Fig. 2. — Numerical (points) and exact (solid lines) solutions att = 20.0. The number of cells is 
m = 128 in the calculations. 
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Fig. 3. — Average deviation from exact solution at t = 20.0 as a function of cell number m. We 
also show reference lines D oc m~^andm~'^ as dashed lines. 
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number density • n [ cm"^ ] 

Fig. 4. — Thermal-equilibrium state of the net cooling function {solid line). Contours of the 
cooling and heating timescales, which correspond to 5.0, 1.0, 0.1, and 0.01 Myr, are shown as 
broken lines. Dotted region shows thermally unstable region, in which the Balbus criterion is 
satisfied. Red and blue lines are the evolutionary tracks of the highest density gas. Red line shows 
the unmagnetized case, which reaches the stable CNM phase at t = 0.5 Myr. Blue line shows the 
magnetized case, which reaches the stable CNM phase at t = 1.6 Myr. 
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Fig. 5. — Density distribution of unmagnetized (to;? Ze/O and magnetized ng/zO cases. Density 
distributions are depicted with the same color scale {bottom). 
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Fig. 6. — Density distribution (top) and magnetic field lines (bottom) of the whole simulation 
domain at the same time depicted in the right panel of Figure [51 Density distributions are depicted 
with the same color scale to Figure [51 Note that the y-scale is expanded in these figures. 



